!/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!*  Copyleft (!c) 2004-2005 Div. Nucl. Med., Cyric, Tohoku Univ.    *
!*  Copyleft (!c) 2005-2011 Elseif Laboratory.                      *
!*  elseifkk@users.sf.net                                           *
!*                                                                  *
!*  All Rights Reversed.                                            *
!*                                                                  *
!*                                                                  *
!*  LICENSE AGREEMENT:                                              *
!*                                                                  *
!*  You agree that we don't agree to anything.                      *
!*                                                                  *
!*                                                                  *
!*  DISCLAIMER OF WARRANTIES:                                       *
!*                                                                  *
!*  This software is provided to you "AS IS," and the authors and   *
!*  contributors disclaim any warranty or liability obligations to  *
!*  you of any kind, implied or express, including but not limited  *
!*  to the warranties of merchantability, fitness for a particular  *
!*  purpose, and non-infringement.                                  *
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
module puty
  implicit none
  type ellipse
     real*8 x(3)
     real*8 a(3)
     real*8 t(3)
     real*8 pmin,pmax
     real*8 dr,r0,rmin
  end type ellipse
  type(ellipse)::e_outer=ellipse(0,(/100.,100.,120./),0,50,100,0.06,0.7,0)
  type(ellipse)::e_inner=ellipse((/0.,-5.,40./),(/70.,75.,60./),(/0.3,0.,0./),150,400,2.5,0.1,0)
  type(ellipse)::e_reb=ellipse((/-35,75,20/),(/20,20,20/),0,10,50,0.04,0.8,0)
  type(ellipse)::e_leb=ellipse((/+35,75,20/),(/20,20,20/),0,10,50,0.04,0.8,0)
  !
  type(ellipse)::f_outer=ellipse(0,(/100.,100.,120./),0,50,150,0.005,0.95,0)
  type(ellipse)::f_outer2=ellipse(0,(/100.,100.,120./),0,200,300,0.01,0.95,0.98)
  type(ellipse)::f_inner=ellipse((/0.,-5.,40./),(/70.,75.,60./),(/0.3,0.,0./),75,200,0.005,0.95,0)
  type(ellipse)::f_reb=ellipse((/-35,75,20/),(/20,20,20/),0,10,100,0.01,0.9,0)
  type(ellipse)::f_leb=ellipse((/+35,75,20/),(/20,20,20/),0,10,100,0.01,0.9,0)
contains
  function rot_x(x,t)
    real*8 rot_x(3)
    real*8,intent(in)::x(3)
    real*8,intent(in)::t
    real*8 c,s
    c=cos(t)
    s=sin(t)
    rot_x(1)=x(1)
    rot_x(2)= c*x(2)+s*x(3)
    rot_x(3)=-s*x(2)+c*x(3)
  end function rot_x
  !
  function rot_y(x,t)
    real*8 rot_y(3)
    real*8,intent(in)::x(3)
    real*8,intent(in)::t
    real*8 c,s
    c=cos(t)
    s=sin(t)
    rot_y(2)=x(2)
    rot_y(3)= c*x(3)+s*x(1)
    rot_y(1)=-s*x(3)+c*x(1)
  end function rot_y
  !
  function rot_z(x,t)
    real*8 rot_z(3)
    real*8,intent(in)::x(3)
    real*8,intent(in)::t
    real*8 c,s
    c=cos(t)
    s=sin(t)
    rot_z(3)=x(3)
    rot_z(1)= c*x(1)+s*x(2)
    rot_z(2)=-s*x(1)+c*x(2)
  end function rot_z
  !
  function rot(x,t)
    real*8 rot(3)
    real*8,intent(in)::x(3)
    real*8,intent(in)::t(3)
    real*8 xx(3)
    xx=rot_x(x,t(1))
    xx=rot_y(xx,t(2))
    rot=rot_z(xx,t(3))
  end function rot
  !
  real*8 function int_norm(e,r)
    type(ellipse),intent(in)::e
    real*8,intent(in)::r
    int_norm=1
  end function int_norm
  !
  real*8 function int_up(e,r)
    type(ellipse),intent(in)::e
    real*8,intent(in)::r
    int_up=1/(1+exp((r-e%r0)/e%dr))
  end function int_up
  !
  real*8 function int_down(e,r)
    type(ellipse),intent(in)::e
    real*8,intent(in)::r
    int_down=1/(1+exp(-(r-e%r0)/e%dr))
  end function int_down
  !
  real*8 function int_pow(e,r)
    type(ellipse),intent(in)::e
    real*8,intent(in)::r
    int_pow=r**e%dr
  end function int_pow
  !
  subroutine set_ellipse(id,a,e,int_f)
    use buty
    type(image_data),intent(inout)::id
    real*8,intent(in)::a(3,4)
    type(ellipse),intent(in)::e
    real*8,external::int_f
    integer*4 ix,iy,iz
    real*8 x(3),r,p,xx(3)
    do iz=1,id%n(3)
       x(3)=iz*id%d(3)-id%offset(3)
       do iy=1,id%n(2)
          x(2)=iy*id%d(2)-id%offset(2)
          do ix=1,id%n(1)
             x(1)=ix*id%d(1)-id%offset(1)
             xx=a.by.x ! global transform
             xx=rot(xx-e%x,e%t)
             r=sqrt(sum((xx/e%a)**2))
             if(r<1.and.r>e%rmin) then
                p=(e%pmax-e%pmin)*int_f(e,r)+e%pmin
                id%img_r8(ix,iy,iz)=p
             end if
          end do
       end do
    end do
    !
  end subroutine set_ellipse
  !
  subroutine mul_random(id,int,flag)
    use buty
    type(image_data),intent(inout)::id
    real*8,intent(in)::int
    integer*4,intent(in)::flag
    real*8,allocatable::r(:,:,:,:)
    allocate(r(id%n(1),id%n(2),id%n(3),2))
    call random_number(r)
    r(:,:,:,1)=sqrt(-2.d0*log(1.d0-r(:,:,:,1)))*sin(2.d0*pi*r(:,:,:,2))
    if(flag==0) then
       id%img_r8=id%img_r8+int*(id%img_r8)*r(:,:,:,1)
    else
       id%img_r8=id%img_r8+sqrt(id%img_r8)*r(:,:,:,1)
    endif
    deallocate(r)    
  end subroutine mul_random
  !
  subroutine emi_head(id,rtp)
    use buty
    type(image_data),intent(inout)::id
    type(rigid_transform_param),intent(in)::rtp
    real*8 a(3,4)
    call set_affine_rigid(rtp,a)
    call set_ellipse(id,a,e_outer,int_up)
    call set_ellipse(id,a,e_inner,int_pow)
    call set_ellipse(id,a,e_reb,int_down)
    call set_ellipse(id,a,e_leb,int_down)
    call mul_random(id,0.2d0,0)
  end subroutine emi_head
  !
  subroutine mri_head(id,rtp)
    use buty
    type(image_data),intent(inout)::id
    type(rigid_transform_param),intent(in)::rtp
    real*8 a(3,4)
    call set_affine_rigid(rtp,a)
    call set_ellipse(id,a,f_outer,int_up)
    call set_ellipse(id,a,f_inner,int_up)
    call set_ellipse(id,a,f_reb,int_up)
    call set_ellipse(id,a,f_leb,int_up)
    call set_ellipse(id,a,f_outer2,int_norm)
    call mul_random(id,0.2d0,1)
  end subroutine mri_head
  !
end module puty
